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We investigate interaction networks that we derive from multivariate time series with methods 
frequently employed in diverse scientific fields such as biology, quantitative finance, physics, earth 
and climate sciences, and the neurosciences. Mimicking experimental situations, we generate time 
series with finite length and varying frequency content but from independent stochastic processes. 
Using the correlation coefficient and the maximum cross-correlation, we estimate interdependencies 
between these time series. With clustering coefficient and average shortest path length, we observe 
unweighted interaction networks, derived via thresholding the values of interdependence, to pos- 
sess non-trivial topologies as compared to Erdos-Renyi networks, which would indicate small-world 
characteristics. These topologies reflect the mostly unavoidable finiteness of the data, which limits 
the reliability of typically used estimators of signal interdependence. We propose random networks 
that are tailored to the way interaction networks are derived from empirical data. Through an 
exemplary investigation of multichannel electroencephalographic recordings of epileptic seizures — 
known for their complex spatial and temporal dynamics — we show that such random networks help 
to distinguish network properties of interdependence structures related to seizure dynamics from 
those spuriously induced by the applied methods of analysis. 

PACS numbers: 



Introduction 

The last years have seen an extraordinary success of 
network theory and its applications in diverse disciplines, 
ranging from sociology, biology, earth and climate sci- 
ences, quantitative finance, to physics and the neuro- 
sciences [1-4]. There is now growing evidence that re- 
search into the dynamics of complex systems profits from 
a network perspective. Within this framework, complex 
systems are considered to be composed of interacting sub- 
systems. This view has been adopted in a large number 
of modeling studies and empirical studies. It is usually 
assumed that the complex system under study can be 
described by an interaction network, whose nodes repre- 
sent subsystems and whose links represent interactions 
between them. Interaction networks derived from empir- 
ical data (multivariate time series) have been repeatedly 
studied in climate science (climate networks, see [5-9] 
and references therein), in seismology (earthquake net- 
works, see, e.g., [10-13]), in quantitative finance (finan- 
cial networks, see e.g. [14-18] and references therein), 
and in the neurosciences (brain functional networks, see 
[19, 20] for an overview). Many interaction networks 
have been reported to possess non-trivial properties such 
as small- world architectures, community structures, or 
hubs (nodes with high centrality), all of which have been 
considered to be characteristics of the dynamics of the 
complex system. 
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When analyzing empirical data one is faced with the 
challenge of defining nodes and inferring links from mul- 
tivariate noisy time series with only a limited number 
of data points due to stationarity requirements. Differ- 
ent approaches varying to some degree across disciplines 
have been proposed. For most approaches, each single 
time series is associated with a node and inference of 
links is based on time series analysis techniques. Bivari- 
ate time series analysis methods, such as the correlation 
coefficient, are used as estimators of signal interdepen- 
dence which is assumed to be indicative of an interaction 
between different subsystems. Inferring links from esti- 
mates of signal interdependence can be achieved in differ- 
ent ways. Weighted interaction networks can be derived 
by considering estimated values of signal interdependence 
(sometimes mapped via some function) as link weights. 
Since methods characterizing unweighted networks are 
well-established and readily available, such networks are 
more frequently derived from empirical data. Besides ap- 
proaches based on constructing minimum spanning trees 
(see, e.g., reference [14]), on significance testing [21-23], 
or on rank-ordered network growth (see, e.g., reference 
[15]), a common practice pursued in many disciplines is 
to choose a threshold above which an estimated value of 
signal interdependence is converted into a link ( "thresh- 
olding", see, e.g., references [5, 12, 16, 20]). Follow- 
ing this approach, the resulting unweighted interaction 
networks have been repeatedly investigated employing 
various networks characteristics, among which we men- 
tion the widely-used clustering coefficient C and average 
shortest path length L to assess a potential small-world 
characteristic, and the node degrees in order to identify 
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hubs. 

As studies employing the network approach grow in 
numbers, the question arises as to how informative re- 
ported results are with respect to the investigated dy- 
namical systems. To address this issue, properties of 
interaction networks are typically compared to those ob- 
tained from network null models. Most frequently, Erdos- 
Renyi random networks [24] or random networks with a 
predefined degree distribution [25, 26] serve as null mod- 
els; network properties that deviate from those obtained 
from the null model arc considered to be characteristic of 
the investigated dynamical system. Only in a few recent 
studies, results obtained from network analyses have been 
questioned in relation to various assumptions underlying 
the network analysis approach. Problems pointed out in- 
clude: incomplete data sets and observational errors in 
animal social network studies [27]; representation issues 
and questionable use of statistics in biological networks 
(see [28] and references therein); challenging node and 
link identification in the neurosciences [29-31]; the issue 
of spatial sampling of complex systems [31-33]. This calls 
not only for a careful interpretation of results but also for 
the development of appropriate null models that incorpo- 
rate knowledge about the way networks are derived from 
empirical data. 

We study - from the perspective of field data analy- 
sis - a fundamental assumption underlying the network 
approach, namely that the multivariate time series are 
obtained from interacting dynamical processes and are 
thus well represented by a model of mutual relationships 
(i.e., an interaction network). Visual inspection of em- 
pirical time series typically reveals a perplexing variety of 
characteristics ranging from fluctuations on different time 
scales to quasi-periodicity suggestive of different types of 
dynamics. Moreover, empirical time series are inevitably 
noisy and finite leading to a limited reliability of esti- 
mators of signal intcrdcpcndencics. This is aggravated 
with the advent of time-resolved network analyses where 
a good temporal resolution often comes at the cost of di- 
minished statistics. Taken together, it is not surprising 
that the suitability of the network approach is notoriously 
difficult to judge prior to analysis. 

We here employ the above-mentioned thresholding- 
approach to construct interaction networks for which we 
estimate signal interdependence with the frequently used 
correlation coefficient and the maximum cross correla- 
tion. We derive these networks, however, from multi- 
variate time series of finite length that are generated by 
independent (non-interacting) processes which would a 
priori not advocate the notion of a representation by 
a model of mutual relationships. In simulation studies 
we investigate often used network properties (clustering 
coefficient, average shortest path length, number of con- 
nected components). We observe that network properties 
can deviate pronouncedly from those observed in Erdos- 
Rcnyi networks depending on the length and the spec- 
tral content of the multivariate time scries. We address 
the question whether similar dependencies can also be 



observed in empirical data by investigating multichan- 
nel electroencephalographic (EEG) recordings of epilep- 
tic seizures that are known for their complex spatial and 
temporal dynamics. Finally, we propose random net- 
works that are tailored to the way interaction networks 
are derived from multivariate empirical time series. 



Methods 

Interaction networks arc typically derived from N mul- 
tivariate time series x% (i G {l,...,iV}) in two steps. 
First, by employing some bivariate time series analysis 
method, interdependence between two time series Xi(t) 
and Xj(t) (t G {1, ...,T}) is estimated as an indicator 
for the strength of interaction between the underlying 
systems. A multitude of estimators [34-40], which differ 
in concepts, robustness (e.g., against noise contamina- 
tions), and statistical efficiency (i.e., the amount of data 
required), is available. Studies that aim at deriving in- 
teraction networks from field data frequently employ the 
absolute value of the linear correlation coefficient to es- 
timate interdependence between two time series. The 
entries of the correlation matrix p c then read 



t=i 

=: |corr(xi,x,-)| , 



(1) 



where x% and a% denote mean value and the estimated 
standard deviation of time series X{ . Another well estab- 
lished method to characterize interdependencies is the 
cross correlation function. Here we use the maximum 
value of the absolute cross correlation between two time 
series, 



max ■ 



with 

£(x<,Xj)(t) 



\/£(zi,Zi)(0)£(x j ,x. ) -)(0) 



£(Xj,Xi)(-T) ,T<0 



(2) 



(3) 



to define the entries of the cross correlation matrix p m . 
As practiced in field data analysis, we normalize the time 
series to zero mean before pursuing subsequent steps of 
analysis. Note that p^ is then the maximum value of 
the absolute cross covariance function. Both interdepen- 
dence estimators are symmetric (p?- = p^ and = p 1 -]) 
and are confined to the interval [0,1]. High values indi- 
cate strongly interdependent time series while dissimilar 
time series result in values close to zero for T sufficiently 
large. 

Second, the adjacency matrix A representing an un- 
weighted undirected interaction network is derived from 
p c (or p m ) by thresholding. For a threshold 6 G [0,1] 
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entries Aij and Aji of A are set to 1 (representing an 
undirected link between nodes i and j) for all entries 
pfj > 9 (pfj > 0, respectively) with i ^ j, and to zero (no 
link) otherwise. In many studies 9 is not chosen directly 
but determined such that the derived network possesses a 
previously specified mean degree k := iV" 1 ki, where 
ki denotes the degree of i, i.e., the number of links con- 
nected to node i. More frequently, 9 is chosen such that 
the network possesses a previously specified link density 
e = k(N — l) -1 . We will follow the latter approach and 
derive networks for fixed values of e. 

To characterize a network as defined by A, a plethora 
of methods have been developed. Among them, the clus- 
tering coefficient C and the average shortest path length 
L are frequently used in field studies. The local clustering 
coefficient C; is defined as 



where 



Ci 



ki(ki-l) ^lj,mAijAj m A mi , if ki > 1 . . 

0, if h e {0,1}. 1 ; 



Ci represents the fraction of the number of existing links 
between neighbors of node i among all possible links be- 
tween these neighbors [1, 2, 41]. The clustering coeffi- 
cient C of the network is defined as the mean of the local 
clustering coefficients, 



C 



1 N 



(5) 



s '■= {(hj) I hj < oo; i,j = 1,...,N} 



(8) 



denotes the set of all pairs (i, j) of nodes with finite short- 
est path. The number of such pairs is given by \S\. Note 
that L ->• for N c -> N. 

In field studies, values of C and L obtained for inter- 
action networks are typically compared with average val- 
ues obtained from an ensemble of random Erdos-Renyi 
(ER) networks [24]. Between every pair of nodes is a 
link with probability e, and links for different pairs exist 
independently from each other. The expectation value 
of the clustering coefficient of ER networks is Cer = e 
[2]. The dependence of the average shortest path length 
Ler of ER networks on e and N is more complicated 
(see references [2, 44]). Almost all ER networks are con- 
nected, if e ^ 1nN/(N — 1). ER networks with a prede- 
fined number of links (and thus link density) can also be 
generated by successively adding links between randomly 
chosen pairs of nodes until the predefined number of links 
is reached. During this process, multiple links between 
nodes arc avoided. 



Results 



Simulation studies 



C quantifies the local interconnectedness of the network 
and a,C e [0,1]. 

The average shortest path length is defined as the av- 
erage shortest distance between any two nodes, 



L := 



N(N ■ 



(6) 



and characterizes the overall connectedness of the net- 
work, kj denotes the length of the shortest path between 
nodes i and j. The definition of the average shortest path 
length varies across the literature. Like some authors, we 
here include the distance from each node to itself in the 
average (la = OVi). Exclusion will, however, just change 
the value by a constant factor of (N + 1)/ (N — 1). 

If a network disintegrates into a number N c of differ- 
ent connected components, there will be pairs of nodes 
(i,j), for which no connecting path exists, in which case 
one usually sets Uj = oo and thus L = oo. In order to 
avoid this situation, in some studies lij in equation (6) is 
replaced by . The quantity defined this way is called 
efficiency [42, 43]. Another approach, which we will fol- 
low here and which is frequently used in field studies, is 
to exclude infinite values of Uj from the average. The 
average shortest path length then reads 



1 1 (i,j)es 



(7) 



We consider time series z,, i £ {1, . . . , N}, whose en- 
tries Zi{t) are drawn independently from the uniform 
probability distribution U on the interval (0, 1). We will 
later study the impact of different lengths T of these ran- 
dom time series on network properties. In order to enable 
us to study the effects of different spectral contents on 
network properties, we add the possibility to low-pass 
filter Zi by considering 



t+M-l 



Xi,M,T 



(t) := M _1 *(0. Zt(l)~U, (9) 



i=t 



where t e {1,...,T}, and 1 < M < T. By definition 
x i i,t(*) = z i{t)^t- With the size M of the moving aver- 
age we control the spectral contents of time series. Wc 
here chose this ansatz for the sake of simplicity, for its 
mathematical treatability, and because the random time 
series with different spectral contents produced this way 
show all properties we want to illustrate. 

In the following we will study the influence of the 
length T of time series on network properties by con- 
sidering Xi_\ t T for different T. For a chosen value of T 
we determine R realizations of x, i t an d we denote each 



realization r with 



,1.T' 



When studying the influence 



of the spectral content we will consider x^m.t' with dif- 
ferent M and with T = 500. Wc chose this value of T' 
because we are interested in investigating time series of 
short length as typically considered in field studies. For a 
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chosen value of M we determine R realizations of Xi.M.T' 

(r) 

and we denote realization r with x\ ^ T , . 

In order to keep the line of reasoning short and clear, 
we will present supporting and more rigorous mathemat- 
ical results in Appendix A and refer to them in places 
where needed. In addition, since we observed most simu- 
lation studies based on p m to yield qualitatively the same 
results as those based on p c , we will present results based 
on p c only and report results of our studies based on p m 
whenever we observed qualitative differences. 

Clustering coefficient 

Let Pij \ t := P c ( x i l,Ti x j\i,t) denote the absolute 
value of the empirical correlation coefficient estimated 

(r) (r) 

for time series x\ { T and x^ [ T , and let us consider R re- 
alizations, r G {1, . . . ,R}- Because of the independence 
of processes generating the time series and because of the 
symmetry of the correlation coefficient, we expect the R 
values of the empirical correlation coefficient calculated 
for finite and fixed T to be distributed around the mean 
value 0. The variance of this distribution will be higher 
the lower we choose T . If we sample one value p^\ T out 

(r) 

of the R values it is almost surely that p\^\ T > 0. Thus 

(r) 

there are thresholds < 9 < p\- \ T for which we would 
establish a link between the corresponding nodes i and j 
when deriving a network. Let us now consider a network 
of N nodes whose links are derived from N time series 
as before. For some 8 > the network will possess links 
and e > 0. We expect to observe e for some fixed 8 > 
to be higher the larger the variance of the distribution of 

(r) 

Pij i t- Likewise, for fixed values of e we expect to find 8 
to be higher the lower we choose a value of T. 

As a first check of this intuition we derive an approx- 
imation e a i for the edge density by making use of the 
asymptotic limit (T — > oo, see Appendix A, Lemma 2 for 
details), 

e a i(0,T)=2$(-\/T0), (10) 

where $ denotes the cumulative distribution function of 
a standard normal distribution. In figure 1 (top left) 
we show the dependence of e a \(8,T) on 8 for exemplary 
values of T. Indeed, e a i((9,T) is decreasing in 8 and T. 

The concession of taking the asymptotic limit when 
deriving equation (10) may limit its validity in the case 
of small values of T in which we are especially interested. 
Thus, we approach this case by simulation studies. Let us 
consider R = 10 6 values of p^J M T obtained for R realiza- 
tions of two time series Xi t M,T, i G {1, 2}, r G {1, . . . , R}. 
We estimate the edge density i(8, M, T) by 

e(6,M,T) := R^^H^^O), (11) 

r 

where H\^ MT {8) = 1 for p^ MT > 8, and else. Note 
that e(8, M, T) does not depend on N. This is because 



e((9, M, T) represents the (numerically determined) prob- 
ability that there is a link between two vertices. The 
dependence of e(8,l,T) on 8 for different values of T 
shown in figure 1 (top left) indicates a good agreement 
between e a i(#, T) and e(8, 1, T) for larger values of T but 
an increasing difference for T < 30. 

We proceed by estimating the clustering coefficient C 
for our model using R realizations of three time series 
Xi,M t T, i G {1, . . . , 3} by 

C(8 M T) "^12, Af,T(^)-^13,M,r(^) -^23, M,t(^) 

(12) 

The dependence of C(8, l,T) on 8 for various T is shown 
in the top right part of figure 1. For fixed T, C(8, 1,T) 
decreases from 1 with increasing values of 8 which one 
might expect due to the decrease of e. However, we also 
observe for 8 > that C(8, 1,T) takes on higher values 
the lower T. 

In order to investigate whether the cluster- 
ing coefficients of our networks differ from those 
of Erdos-Rcnyi networks we use equation (11) 
and obtain CW(e) = C(6{e, M,T), M,T) with 
8(e,M,T) = inf{<9 : e(8,M,T) > e}. This allows 
the comparison with Cer(c) = £ by considering the 
ratio 7(e, M,T) := Cm,t{^)/Cer(^)- Remarkably, 
7(e, 1, T) ^> 1 for a range of values of e and T (cf. lower 
left part of figure 1). 7(e, 1,T) even increases for small e 
and T. This indicates that there is a relevant dependence 
between the three random variables p%j,M,T^ Pu,m,t, 
and Pji t M,T for different indices I and small T. For 
T — > oo and fixed edge density, C converges to Cer, 
because the dependence between the random variables 
Pij,M,T, i,j G {l,...,iV}, vanishes (i.e., the random 
variables will converge in distribution to independent 
normal random variables). 

In order to gain deeper insights into the influence 
of the spectral contents of random time series on the 
clustering coefficient, we repeat the steps of analysis 
with time series Xi t M,T> , where T" = 500 is kept fix, 
and we choose different values of M. Figure 1 (top 
panels and lower left) shows that the higher the amount 
of low-frequency contributions (large M) the higher 
e(6,M,T') and 6(6, M,T') (for 8 > 0), and the higher 
j(e,M,T') (for e < 1). The difference between Erdos- 
Renyi networks and our time series networks becomes 
more pronounced (j(e, M,T') 1) the smaller e and 
the higher M. 

Given the similar dependence of 7, C, and e on T and M, 
we hypothesize that the similarity can be traced back to 
similar variances of pij.\,T and Pij,M,T' for some values of 
T and M. By making use of the asymptotic variance of 
the limit distributions of T — > 00, we derive an expression 
relating Yax(pij \jr) an d Var(/3y Af.r) to each other (see 
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FIG. 1: Simulation results for the edge density, the clustering coefficient, and the effective length. Top row: De- 
pendence of edge density e(9, M, T) (left) and of clustering coefficient C{9, M, T) (right) on the threshold 6 for different values of 
the size M of the moving average and of the length T of time series. Values of edge density e a i(#, T) obtained with the asymptotic 
limit (equation (10)) are shown as lines (top left). Bottom left: Dependence of the ratio -y(e,M,T) = Cm,t(£)/Cer.(£) on edge 
density e. Note, that we omitted values of estimated quantities obtained for 6 £ {9 : (RT 1 t(^)^i3 m t(^)) < 1CP 3 } 

since the accuracy of the statistics is no longer guaranteed. Bottom right: Dependence of effective length T e g as determined 
by equation (14) (black line) and its numerical estimate T e ff (red markers) on M. 



Appendix A, Lemma 1), 

Var(pij j Af,r) ~ ff(M)Vax(py,i,r), 

with<7(M) = ^M + ^, (13) 

We are now able to define an effective length T c s of time 
series, 



T cff (M) 



T 



(14) 



for which Var(py i i ) T ttB ) ~ Var( ( Oy.M.T')- I n the lower 
right part of figure 1 we show T c g (M) in dependence on 
M. To investigate whether equation (14) also holds for 
small values of T, we determine numerically, for differ- 
ent values of 9, C(0, 1,T) for T e {3,...,T'} as well 



as C(0, M, T') for some chosen values of M. Eventu- 
ally, we determine for each value of M a value of T, for 
which 0(6,1, T) and C(9, M, T") curves match in a least- 
squares sense, and denote this value with T e g (see the 
lower right part of figure 1). We observe a maximum de- 



viation 



Toff — T e ff 



2 and conclude that equation (14) 

indeed holds for small length T of time series. Moreover, 
numerically determined dependencies of £ on fl, C on f?, 
as well as 7 on e for pairs of values (M, T') show a re- 
markable similarity to those dependencies obtained for 
pairs of values (1, T c g). 

Thus, the clustering coefficient of networks derived 
from random time series of finite length and/or with 
a large amount of low-frequency contributions is higher 
than the one of Erdos-Rcnyi (ER) networks - indepen- 
dently of the network size N (cf. equation (12)). This 
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FIG. 2: Simulation results for the average shortest path length. Dependence of the average shortest path length 
L(e,M,T) (left) and of the ratio A(e, M, T) = L(e, 1, T)/Ler(e) (right) on edge density e for different values of the size M of 
the moving average and of the length T of time series. Lines are for eye-guidance only. 



difference becomes more pronounced the lower the edge 
density e, the lower the length T of time series, and the 
larger the amount of low-frequency contributions. These 
results point us to an important difference between ER 
networks and our model networks: possible edges in ER 
networks are not only (1) equally likely but also (2) inde- 
pendently chosen to become edges. While property (1) 
is fulfilled in our model networks, property (2) is not. 



Average shortest path length 

Next we study the impact of the length of time series 
and of the amount of low-frequency contributions on the 
average shortest path length L of our model networks 
by employing a similar but different simulation approach 
as used in the previous section. To estimate L, we con- 
sider R = 100 networks with a fixed number of nodes 
(JV = 100). We derive our model networks by thresh- 
olding pV MtT , i,j € {1,...,N}, r G {1,...,R}. Let 

L,( r '(e,l,T) denote the average shortest path length for 
network r with M = 1 and different values of T, and 
L {r \e,M,T') the average shortest path length for net- 
work r with fixed value of T (T = T' = 500) and different 
values of M. With L^(e) we refer to the average short- 
est path length obtained for the r-th ER network of size 
N and edge density e. Mean values over realizations will 
be denoted as L(e, 1, T), L(e, M, T"), and Ler(c) respec- 
tively. Finally, we define A(e,l,T) := L(e, l,T)/L ER (e) 
and A(e, M, T') := L(e, M, T')/L ER (e). 

In figure 2 we show the dependence of L and A on e for 
various values of M and T. All quantities decrease as e 
increases which can be expected due to additional edges 
reducing the average distances between pairs of nodes of 
the networks. For fixed e<l, L takes on higher values 



the higher M or the lower T. With equation (14) we have 
L(e, 1, T c ff) ps L(e, A/, T") which resembles the results ob- 
tained for the clustering coefficient. Differences between 
the average shortest path lengths of our model networks 
and ER networks (as characterized by A) become more 
pronounced the higher M and the lower T. For edge den- 
sities typically reported in field studies (e ss 0.1), how- 
ever, differences are less pronounced (A < 1.2, cf. figure 2 
right) than the ones observed for the clustering coefficient 
(7 > 2 for selected values of M and T, cf. figure 1 bot- 
tom left). We obtained qualitatively similar results for 
small (N = 50) and large numbers of nodes (N = 500). 

Number of connected components and degree distribution 

Since the number of connected components of a given 
network might affect network characteristics such as the 
average shortest path length (see equation (7)), we in- 
vestigate the impact of different length of time series and 
of the amount of low-frequency contributions on the av- 
erage number of connected components N c (e, M,T) of 
the networks derived from x^i^tt an< ^ x i.M.T'- We de- 
termine N c (e, M,T) as the mean of R realizations of 
the corresponding networks and with iV Cj ER(e) we de- 
note the mean value of the number of connected com- 
ponents in R realization of ER networks. For the edge 
densities considered here we observe ER networks to be 
connected (cf. figure 3), -/V C; er ~ 1, which is in agree- 
ment with the connectivity condition for ER networks, 
e > bxN/(N - 1) ps 0.05 (for N = 100). Similarly, 
we observe N c (e,l,T e g) ps 1, even for small values of 
T (cf. figure 3 right). In contrast, N c {e,M,T') takes 
on higher values the lower e and the higher M (cf. fig- 
ure 3 left). In order to achieve a better understanding 
of these findings, we determine degree probability distri- 
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FIG. 3: Simulation results for the number of connected components. Dependence of the number of connected 
components N c (e, M,T) on the edge density e for different values of the size M of the moving average (left, for T = 500) and 
of the length T of time series (right, for M = 1). Lines are for eye-guidance only. 



butions of our model networks. Let pk denote the esti- 
mated probability of a node to possess a degree k, i.e., 
Pk = #{i {r) : 4 r) = k,r S {1,...,R}}/(NR). With 
Pfc(e, M, T) we will denote the estimated degree distribu- 
tion for networks derived from x^m,t- We briefly recall 
that the degree distribution of ER networks pfc^ER fol- 
lows a binomial distribution, 

which we show in figure 4 for N = 100 and various e 
(top panels and lower left panel). In the same figure we 
present our findings for pk(e, M,T) for various values of 
T = T Q g and M. We observe Pfe(e, 1, T c g) to be equal to 
_Pfc,_/v,ER.(e) within the error to be expected due to the lim- 
ited sample size used for the estimation. For f>k(e, M, T"), 
however, we observe striking differences in comparison to 
the previous degree distributions. In particular, for de- 
creasing e and higher M, the probability of nodes with 
degree k = increases, which leads to networks with dis- 
connected single nodes, thereby increasing the number of 
connected components of the network. 

We hypothesize that the observed differences in the 
number of connected components as well as in the de- 
gree distributions are related to differences in the spec- 
tral content of different realizations of time series x[ M t> 
for M > 1. In particular, a node i with a low degree 

(r) 

ki might be associated with a time series x\ M T , , which 
possesses, by chance, a small relative amount of low fre- 
quency contributions (or, cquivalcntly, a large relative 
amount of high frequency contributions). 

In order to test this hypothesis, we generate R rcaliza- 

(r) 

tions of N = 100 time series x^ ^ j,, and estimate their 

periodograms P^ljif) for frequencies / € {0, . . . , /Nyq} 
using a discrete Fourier transform [45]. /Nyq denotes 



the Nyquist frequency, and periodograms are normalized 

such that T,fPi%(f) = !• From the same time scries, 
we then derive the networks using e = 0.1 and determine 

(r) 

the degrees k\ ' . For some fixed / > we define the to- 
tal power above /' (upper frequency range) as P^m = 

Ep™ ^ ( ,m(/)> and the total P° wer below /' ( lower fre- 
quency range) as P^m = ^2f=o ^ m(/)- For eacn rea h 
ization r we estimate the correlation coefficients between 
the degrees and the corresponding total power contents in 
upper and lower frequency range, = corr(fc^, P^}^) 
and = corr(fc( r ), P^'^), respectively, and deter- 

mine their mean values, Kl(M) = R~ 1 J2r K h' > ano - 

k h (M) = R- 1 J2 r k h ) - Note that «l(M) = -« H (M) by 
construction. We choose /' = f'(M) such that 40% of 
the total power of the filter function associated with the 
moving average is contained within the frequency range 

/e [0,/']. 

For increasing M we observe in the lower right panel 
of figure 4 the degrees to be increasingly correlated with 
Pm ' which corresponds to an anti-correlation of de- 
grees with P^ • Thus, as hypothesized above, the ob- 
served differences in the degree distributions can indeed 
be related to the differences in the power content of the 
time series. We mention that the exact choice of /' does 
not sensitively affect the observed qualitative relation- 
ships as long as < /' <C /Nyq is fulfilled. 

We briefly summarize the results obtained so far, which 
indicate a striking difference between networks derived 
from independent random time series using p c or p m (cf. 
equations (1) and (2)) and corresponding ER networks. 
First, we observed the clustering coefficient C and the 
average shortest path length L of our networks to be 
higher the lower the length T of the time series (cf. fig- 
ures 1 and 2). Second, for some fixed T we observed 
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FIG. 4: Simulation results for the degree distribution, (a-c) Degree distributions pk(e, M,T) estimated for R = 1000 
realizations of networks derived from time series Xi,M,T {N = 100) via thresholding using various edge densities e = k(N — 1) _1 
and for selected values of the size M of the moving average and of the length T of time series. The symbol legend in (a) also 
holds for (b) and (c). (d) Dependence of correlation (kl(M)) between node degrees and power content in the lower frequency 
range on the size M of the moving average. Mean values of correlations obtained for R — 100 realizations of networks for each 
value of M are shown as crosses and standard deviations as error bars. Stars indicate significant differences in comparison to 
(Bonferroni corrected pair-wise Wilcoxon rank sum tests for equal medians, p < 0.01). Lines are for eye-guidance only. 



C and L to be higher the larger the amount of low fre- 
quency components (as parametrized by M) in the time 
series. In addition, these contributions led to an increas- 
ing number of connected components in our networks 
and to degree distributions that differed strongly from 
those of the corresponding ER networks (cf. figures 3 
and 4) . We mention that L as defined here (cf. equation 
(7)) tends to decrease for networks with an increasing 
number N c of connected components, and L — >• for 
N c — > N. L thus depends non-trivially on the amount of 
low frequency components in the time series. Third, for 
small edge densities e and for short time series lengths 
or for a large amount of low frequency components, the 
clustering coefficient deviates more strongly from the one 
of corresponding ER networks (7 > 2) than the average 
shortest path length (A < 1.2; cf. figure 2 right and fig- 
ure 1 (bottom left)). Networks derived from independent 
random time series can thus be classified as small world 
networks if one uses 7> 1 and A ~ 1 as practical crite- 
rion, which is often employed in various field studies (cf. 



[31] and references therein). 



Field data analysis 

The findings obtained in the previous section indicate 
that strong low frequency contributions affect network 
properties C and L in a non-trivial way. We now inves- 
tigate this influence in electroencephalographic (EEG) 
recordings of epileptic seizures that are known for their 
complex spatial and temporal changes in frequency con- 
tent [46-49]. We analyze the multichannel (N — 53 ± 21 
channels) EEGs from 60 patients capturing 100 epilep- 
tic seizures reported in reference [50]. All patients had 
signed informed consent that their clinical data might be 
used and published for research purposes. The study pro- 
tocol had previously been approved by the ethics com- 
mittee of the University of Bonn. During the prcsur- 
gical evaluation of drug-resistant epilepsy, EEG data 
were recorded with chronically implanted strip, grid, or 
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FIG. 5: Evolving relative amount of power during epileptic seizures. (Left) Relative amount of power P contained 
in the S- (Ps, black), (Pa, blue), a- (P a , green), and /3- (Pp, red) frequency bands during an exemplary seizure. Profiles 
are smoothed using a four-point moving average. Grey-shaded area marks the seizure. (Right) Mean values (Ps, P$, P a , Pp) 
of the relative amount of power averaged separately for pre-seizure, discretized seizure, and post-seizure time periods of 100 
epileptic seizures. Lines are for eye-guidance only. 



depth electrodes from the cortex and from within rele- 
vant structures of the brain. The data were sampled at 
200 Hz within the frequency band 0.5 — 70 Hz using a 16- 
bit analog-to-digital converter. Electroencephalographic 
seizure onsets and seizure ends were automatically de- 
tected [51], and EEGs were split into consecutive non- 
overlapping windows of 2.5 s duration (T = 500 sampling 
points). Time series of each window were normalized to 
zero mean and unit variance. We determined p c and p m 
for all combinations of EEG time series from each window 
and derived networks with a fixed edge density e = 0.1 in 
order to exclude possible edge density effects. With L c 
and C c as well as L m and C m we denote characteristics 
of networks based on p c and p m , respectively. In order 
to simplify matters, we omit the window indexing in the 
following. 

We investigate a possible influence of the power content 
of EEG time series on the clustering coefficient and the 
average shortest path length by comparing their values to 
those obtained from ensembles of random networks that 
are based on properties of the EEG time series at two 
different levels of detail. For the first ensemble and for 
each patient we derive networks from random time series 
with a power content that approximately equals the mean 
power content of all EEG time series within a window. 
Let Pi(f) denote the estimated periodogram of each EEG 
time series i, and with P(f) = N^ 1 Pi(f) we denote 
the mean power for each frequency component / over all 
N time series. We normalize P(f) such that P(f) 

. We generate N random time series of length T = 500 
whose entries are independently drawn from a uniform 
probability distribution, and we filter these time series 
in the Fourier domain using P(f) as filter function. 
We normalize the filtered time scries to zero mean and 
unit variance and derive a network based on p c or p m 



(e = 0.1). We use 20 realizations of such networks per 
window in order to determine the mean values of network 
characteristics C^ 1 ' and as well as c£' and 
based on p c or p m , respectively. Since the power spectra 
of all time series equal each other, these random networks 
resemble the ones investigated in the previous section. 

With the second ensemble, we take into account that 
the power content of EEG time series recorded from dif- 
ferent brain regions may differ substantially. For this 
purpose we make use of a well established method for gen- 
erating univariate time series surrogates [52, 53], which 
have power spectral contents and amplitude distributions 
that are practically indistinguishable from those of EEG 
time series but are otherwise random. Amplitudes are 
iteratively permuted while the power spectrum of each 
EEG time series is approximately preserved. Since this 
randomization scheme destroys any significant linear or 
non-linear dependencies between time series, it has been 
successfully applied to test the null hypothesis of inde- 
pendent linear stochastic processes. For each patient, we 
generated 20 surrogate time series for each EEG time 
series from each recording site and each window, and 
derived networks based on cither p c or p m (e = 0.1). 
Mean values of characteristics of these random networks 



arc denoted as and L^' as well as Cm' and L„ 
respectively. 

We begin with an exemplary recording of a seizure of 
which we show in figure 5 (left) the temporal evolution 
of the relative amount of power in the 5- (0-4 Hz, Ps), 
i?- (4-8 Hz, P#), a- (8-12 Hz, P a ), and f3- (12-20 Hz, 
Pp) frequency bands. Prior to the seizure the (5-band 
contains more than 50% of the total power which is then 
shifted towards higher frequencies and back towards low 
frequencies at seizure end. P$ is even higher after the 
seizure than prior to the seizure. 
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FIG. 6: Evolving network properties during an exemplary epileptic seizure. Network properties C c and L c (top row, 
black lines) as well as C m and L m (bottom row, black lines) during an exemplary seizure (cf. figure 5 (left)). Mean values and 

/ 2 ^ f21 (2^ ("2") 

standard deviations of network properties obtained from surrogate time series (Cc , L c , C m , im ) are shown as blue lines 
and blue shaded areas, respectively, and mean values and standard deviations of network properties obtained from the overall 
power content model (Cc 1 ', L c , Cm\ Lm ) are shown as red lines and red shaded areas, respectively. Profiles are smoothed 
using a four-point moving average. Grey-shaded area marks the seizure. For corresponding Erdos-Renyi networks Cer 
and Z/er ~ 2.4 for all time windows. 
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In figure 6 we show the temporal evolution of network 
properties obtained for this recording based on p c (top 
panels) and p m (bottom panels) . During the seizure both 
the clustering coefficients C c and C m and the average 
shortest path lengths L c and L m show pronounced dif- 
ferences to the respective properties obtained from the 
random networks. These differences are less pronounced 
prior to and after the seizure, where C$ and 1^' even 
approach the values of C m and L m , respectively. Cc 

and Cm decrease during the seizure and already increase 
prior to seizure end, resembling the changes of Pg (cf. 
figure 5 (left)). This is in accordance with results of 
our simulation studies: there we observed the cluster- 
ing coefficient to be higher the larger the amount of low 
frequency components in the time series; this could also 
be observed, but to a much lesser extent, for the av- 
erage shortest path length. Indeed, l£' and vary 
little over time, and is only slightly increased after 
the seizure, reflecting the high amount of power in the 



(5-band. 

We only observe small deviations between Cc 1 '' and 
as well as between L^ 1 ' and L c 2 \ which appear 
to be systematic (for many windows Cc <Cc and 
L C 1 ' I >L C 2 ' 1 ). These suggest that for interaction networks 
derived from p c , both random network ensembles appear 
appropriate to characterize the influence of power in low 
frequency bands on clustering coefficient and the average 
shortest path length. In contrast, we observed differ- 
ences between ci 1 ' and Cm , as well as between Lm and 
Lm ■ These differences were most pronounced during the 
seizure and for Lm and Lm also after the seizure. This 
finding indicates that the clustering coefficient and aver- 
age shortest path length of interaction networks derived 
from p m depend sensitively on the power contents of EEG 
time series recorded from different brain regions. Thus, 
for these interaction networks only the random networks 
that account for the complex changes in frequency con- 
tent of different brain regions prior to, during, and after 
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FIG. 7: Evolving network properties averaged over 100 epileptic seizures. Mean values (black) of network properties 
C c (top left), L c (top right), C m (bottom left), and L m (bottom right) averaged separately for pre-seizure, discretized seizure, 
and post-seizure time periods of 100 epileptic seizures. Mean values of corresponding network properties obtained from the first 
and the second ensemble of random networks are shown as red and blue lines, respectively. All error bars indicate standard 
error of the mean. Lines are for eye-guidance only. 



seizures appear appropriate to characterize the influence 
of power in low frequency bands on clustering coefficient 
and the average shortest path length. 

We continue by studying properties of networks de- 
rived from the EEG recordings of all 100 focal onset 
seizures. Due to the different durations of seizures (mean 
seizure duration: 110 ± 60 s) we partitioned each seizure 
into 10 equidistant time bins (see reference [50] for de- 
tails) and assigned the time-dependent network proper- 
ties to the respective time bins. For each seizure we in- 
cluded the same number of pre-seizure and post-seizure 
windows in our analysis and assigned the correspond- 
ing time-dependent network properties to one pre-seizure 
and one post-seizure time bin. Within each time bin we 
determined the mean value (e.g., C c ) and its standard 



error for each property. In figure 5 (right), we show for 
each time bin the mean values of the relative amount of 
power in different frequency bands of all seizure record- 
ings (-P<s, P&, P a , Pp). Similar to the exemplary record- 
ing (cf. figure 5 (left)), we observe a shift in the relative 
amount of power in low frequencies prior to seizures to- 
wards higher frequencies during seizures and back to low 
frequencies at seizure end. The amount of power in the 
(5-band is on average higher in the post-seizure bin than 
in the pre-seizure bin. 

In figure 7 we show the mean values of properties of 
networks in each time bin for all seizures. We observe 
C c (1) , C { c 2 \ L^, Li 2 \ Ci 1 ', and L$ to decrease dur- 
ing seizures and to increase prior to seizure end thereby 
roughly reflecting the amount of power contained in low 
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FIG. 8: Evolving normalized network properties averaged over 100 epileptic seizures. Mean values of Cc/Cq 2 ' and 

Cm/C$ (left) as well as Lc/L^ and Lm/iS (right) averaged separately for pre-seizure, discretized seizure, and post-seizure 
time periods of 100 epileptic seizures. All error bars indicate standard error of the mean. Lines are for eye-guidance only. 



frequencies (cf. figure 5 (right), P$). Ci and Cc Z; and to 

a lesser extent also Z c and L c roughly follow the same 
course in time, however, with a slight shift in the range 
of values as already observed in the exemplary record- 
ing of a seizure (cf. figure 6). Differences between both 
random network ensembles are most pronounced in net- 
work properties based on p m , i.e., between Cm and Cm 
as well as between and Zm . This corroborates the 
observation that the clustering coefficient and the aver- 
age shortest path length of the random networks based 
on p m depend more sensitively on the power contents 
of EEG time series recorded from different brain regions 
than the respective quantities derived from p c . While L c 
and L m show a similar course in time, reaching a maxi- 
mum in the middle of the seizures, we observe a remark- 
able difference between C c and C m prior to end of the 
seizures, where the amount of power in low frequencies is 
large. While C m decreases at the end of the seizures, C c 
does not and remains elevated after seizures. Interest- 
ingly, considering the corresponding quantities obtained 
from the second random network ensemble, C m fluctu- 
ates around 0.3 ± 0.01 and does not increase at the end 
of seizures, while, in contrast, Cc increases at the end 
of the seizures, traversing an interval of values roughly 
three times larger than the interval containing values of 
C m ; . Taken together these findings suggest that the pro- 
nounced changes of the frequency content of EEG time 
series seen during epileptic seizures influence the values 
of the clustering coefficient and the average shortest path 
length. 

A comparison of some value of a network property with 
the one obtained for a random network with the same 
edge density and number of nodes is typically achieved 
by calculating their ratio. If ER networks are used for 
comparison, the value of a network property is rescaled 
by a constant factor. In this case, the time-dependent 
changes of network properties shown in figure 7 will be 
shifted along the ordinate only. In order to take into 



i(2) 



account the varying power content of EEG time series 
recorded from different brain regions we instead normal- 
ize the clustering coefficients and the average shortest 
path lengths with the corresponding quantities from the 
second random network ensemble Cm , and 

L m (cf. figure 8). We observe the normalized network 
properties to describe a concave-like movement over time 
indicating a reconfiguration of networks from more ran- 
dom (before seizures) towards a more regular (during 
seizures) and back towards more random network topolo- 
gies. This is in agreement with previous observations 
using a different and seldom used thresholding method 
[50]. 



Discussion 

The network approach towards the analysis of empir- 
ical multivariate time series is based on the assumption 
that the data is well represented by a model of mutual re- 
lationships (i.e., a network). We studied interaction net- 
works derived from finite time series generated by inde- 
pendent processes that would not advocate a representa- 
tion by a model of mutual relationships. We observed the 
derived interaction networks to show non-trivial network 
topologies. These are induced by the finiteness of data, 
which limits reliability of estimators of signal interdepen- 
dence, together with the use of a frequently employed 
thresholding technique. Since the analysis methodology 
alone can already introduce non-trivial structure in the 
derived networks, the question arises as to how informa- 
tive network analysis results obtained from finite empir- 
ical data are with respect to the studied dynamics. This 
question may be addressed by defining and making use 
of appropriate null models. In the following, we briefly 
discuss two null models that arc frequently employed in 
field studies. 

Erdos-Rcnyi (ER) networks represent one of the ear- 
liest and best studied network models in mathematical 
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literature and can be easily generated. They can be used 
to test whether the network under consideration complies 
with the notion of a random network in which possible 
edges are equally likely and independently chosen to be- 
come edges. We observed that clustering coefficient C 
and average shortest path length L for interaction net- 
works derived from finite random time series differed pro- 
nouncedly from those obtained from corresponding ER 
networks, which would likely lead to a classification of 
interaction networks as small-world networks. Since the 
influence of the analysis methodology is not taken into 
account with ER networks, they may not be well suited 
for serving as null models in studies of interaction net- 
works derived from finite time series. 

Another null model is based on randomizing the net- 
work topology while preserving the degrees of nodes 
[26, 54, 55]. It is used to evaluate whether the network 
under consideration is random under the constraint of a 
given degree sequence. Results of our simulation stud- 
ies point out that the structures induced in the network 
topology by the way how networks are derived from em- 
pirical time series cannot be related to the degree se- 
quence only. We observed that C and L from interaction 
networks remarkably depended on the finitcness of the 
data, while the degree distribution did not (cf. figure 4 
(a-c), M = 1). The usefulness of degree-preserving ran- 
domized networks has also been subject of debate since 
they do not take into account different characteristics of 
the data and its acquisition [56, 57]. Moreover, the link- 
switching algorithm frequently employed for generating 
such networks has been shown to non-uniformly sample 
the space of networks with predefined degree sequence 
(see, e.g., references [25, 58]). This deficiency can be ad- 
dressed by using alternative randomization schemes (see, 
e.g., [58-60] and references therein). 

We propose to take into account the finite length and 
the frequency contents of time series when defining null 
models. For this purpose we applied the same method- 
ological steps as in field data analysis (estimation of sig- 
nal interdependence and thresholding of interdependence 
values to define links) but used surrogate time series [53] 
to derive random networks (second ensemble). These 
surrogate time scries comply with the null hypothesis 
of independent linear stochastic processes and preserve 
length, frequency content, and amplitude distribution of 
the original time series. For these random networks, we 
observed (in our simulation studies) dependencies be- 
tween properties of networks and properties of time se- 
ries: the clustering coefficient C, and, to a lesser extent, 
the average shortest path length L are higher the higher 
the relative amount of low frequency components, the 
shorter the length of time series, and the smaller the edge 
density of the network. Results obtained from an anal- 
ysis of interaction networks derived from multichannel 
EEG recordings of one hundred epileptic seizures con- 
firm that the pronounced changes of the frequency con- 
tent seen during seizures influence the values of C and 
L. Comparing these network characteristics with those 



obtained from our random networks allowed us to distin- 
guish aspects of global network dynamics during seizures 
from those spuriously induced by the applied methods of 
analysis. 

Our random networks will likely be classified as small- 
world networks when compared to ER networks which 
might indicate that small-world topologies in networks 
derived from empirical data as reported in an ever in- 
creasing number of studies can partly or solely be related 
to the finite length and frequency content of time series. 
If so, small-world topologies would be an overly com- 
plicated description of the simple finding of finite time 
series with a large amount of low frequency components. 
In this context, our approach could be of particular in- 
terest for studies that deal with short time series and low 
frequency contents, as, for example, is the case in rest- 
ing state functional magnetic resonance imaging studies 
(see, e.g., references [61-65]). In such studies, taking into 
account potential frequency effects could help to unravel 
information on the network level that would be otherwise 
masked. 

We observed the degrees of nodes of our random net- 
works to be correlated with the relative amount of power 
in low-frequencies in the respective time series (cf. fig- 
ure 4). The degree of a node has been used in field studies 
as an indicator of its centrality in the network (see, e.g., 
[2, 66] and references therein). Particular interest has 
been devoted to nodes which are highly central (hubs). In 
this context it would be interesting to study whether find- 
ings of hubs in interaction networks can partly or solely 
be explained by the various frequency contents of time 
series entering the analysis. In such a case, hubs would be 
a complicated representation of features already present 
on a single time series level. We are confident that our 
random networks can help to clarify this issue. 

Our simulation studies were based on the simplified as- 
sumption that power spectra of all time series from which 
a network is derived are approximately equal. The de- 
pendencies of C and L on the power content could also 
be observed qualitatively for networks derived from EEG 
time series - that were recorded from different brain re- 
gions and whose power spectra may differ substantially 
among each other - but only if link definition was based 
on thresholding the values of the correlation coefficient 
(p c ). Thus, estimating mean power spectra of multivari- 
ate time series can provide the experimentalist with a 
rule of thumb for the potential relative increase of C and 
L in different networks based on the correlation coeffi- 
cient. This rule of thumb, however, might not be helpful 
if the maximum value of the absolute cross correlation 
(p m ) is used to estimate signal interdependencies. In this 
case, C and L depended sensitively on the heterogeneity 
of power spectra (see the second random network en- 
semble). It would be interesting to investigate in future 
studies, which particular properties of p° and p m can be 
accounted for these differences. 

We close the discussion with two remarks, the first 
being of interest for experimentalists. Our findings also 
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shed light on a network construction technique that relies 
on significance testing in order to decide upon defining 
a link or not [21]. For this purpose, a null distribution 
of a chosen estimator of signal interdependence (p m ) is 
generated for each pair of time series and a link is es- 
tablished if the null hypothesis of independent processes 
generating the time series can be rejected at a predefined 
significance level. It was suggested in Ref. [21] to use a 
limited subset of time series in order to minimize compu- 
tational burden when generating null distributions. Our 
findings indicate that networks constructed this way will 
yield an artificially increased number of false positive or 
of false negative links which will depend on the frequency 
contents of time series being part or not part of the sub- 
set. Our second remark is related to network modeling. 
By choosing some threshold and generating time series 
that satisfy the relation between the size of the moving 
average and the length of time series, networks can be 
generated which differ in their degree distributions but 
approximately equal in their clustering coefficient and 
average shortest path length. This property could be of 
value for future modeling studies. 

To summarize, we have demonstrated that interaction 
networks, derived from finite time series via threshold- 
ing an estimate of signal interdependence, can exhibit 
non-trivial properties that solely reflect the mostly un- 
avoidable finiteness of empirical data, which limits the 
reliability of signal interdependence estimators. Address- 
ing these influences, we proposed random network models 
that take into account the way interaction networks are 
derived from the data. With an exemplary time-resolved 
analysis of the clustering coefficient C and the average 
shortest path length L of interaction networks derived 
from multichannel clcctroencephalographic recordings of 
one hundred epileptic seizures, we demonstrated that our 
random networks allow one to gain deeper insights into 
the global network dynamics during seizures. Here we 
concentrated on C and L but we also expect other net- 
work characteristics to be influenced by the methodolo- 
gies used to derive interaction networks from empirical 
data. Analytical investigations of properties of our ran- 
dom networks and the development of formal tests for de- 
viations from these networks may be regarded as promis- 
ing topics for further studies. Other research directions 
are related to the framework we proposed to generate 
random networks from time series. For example, parts of 
the framework may be exchanged in order to study net- 
work construction methodologies other than threshold- 
ing (e.g., based on minimum spanning trees [14] or based 
on allowing weighted links) or other widely used linear 
and nonlinear methods for estimating signal interdepen- 
dence [34, 35, 39]. Other surrogate concepts [67-72] may 
allow for defining different random networks tailored to 
various purposes. We believe that research into network 
inference from time series and into random network mod- 
els that incorporate knowledge about the way networks 
are derived from empirical data can decisively advance 
applied network science. This line of research can con- 
tribute to gain a better understanding of complex dy- 
namical systems studied in various scientific fields. 



Appendix A 
Lemma 1 

For every i,j £ {1, . . . , N} with i ^ j, we have the follow- 
ing limit of the probability distribution of the empirical 
correlation: 



P 



with 



T 



9{M) 



<-"n(.r, ;i; / ..(', ;;/ 7 ; < x ^ $(a;) (Al) 



g(M) 



2 , r 1 1 

3 M+ 3M 



as T — >• oo, where $ denotes the cumulative distribution 
function of a standard normal random variable. 
Proof. In order to simplify the presentation, we write 
Vi,M,T(t) = x. hM ,T(t) - §, so that Ey iiM ,T(t) = 0. First 
note that yi t M,r{t) is a M-dependent sequence, i.e. for 
\s — t\ > M, yi,M,T{s) and yi,M,r{t) are independent. So 
we have that the covariance 



Cov{yi,M,T{^)Vi,M,T{l),yiM,T{t)yj,M,T{t)) = 

for T> M. Additionally, 

Cov (yi,M,T(l)yj,M,T(l),yi,ALT(t)yj,M,T(t)) = 

Cov (yj,M,T(l), Vi,M,T{t)) Cov (j/j,M,T(l)) yj,M,T(t)) 

(A2) 

and Cov (zi(s), Zi(t)) = Var (zi(l)) if s = t and otherwise 
Cov (zi(s),Zi(t)) = 0. For 1 < t < M, we obtain by the 
definition of the moving average and the independence of 
the underlying process Zj(t), t € N that 

COV (yj,Af,T(l)%,M,T(l), yi,M,T(t)Vj,M,T(t)) (A3) 
(M-(t-l) \ 2 

-A £ Var (*,(.))! (A4) 



M 4 



(M - (t - \)Y\bx z (*(1)) . (A5) 



By the central limit theorem for M-dependent random 
variables, see reference [73], 



Var ( i J2t=i yi,M,T(t)yj,M,r(t) 



1 T 

^HwWtoTfi) (A6) 



converges in distribution to a standard normal random 
varibalc as T — > oo. Furthermore, we have the following 
convergence for the variance as T — > oo: 
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M 

-> Var(y i) A ? r,T(l)yj,Ar,T(l)) + 2 E Cov (^,m,t(1)%,m,t(1), Vi,M,T{t)yj,M,T{t)) 



t=2 



(«(!)). (A7) 



The last equality follows easily by E"=i * 2 = an( i T Et=i Vi,M,T(t) — > 0, so we get 



ra(n+i)(2w+i) ^ With the same central limit theorem, 

6 



Et=i Vi,M,T(t) converges to a normal limit, so 
^t- ELi Vi,M,T(t) -> in probability and consequently 1 ^(^.Af.rft) - Z/i.m.t) 



t=i 

v 2 



V t=i /V t=i / t=i V t=i 

= (^H(^H 4 " (A8) ^Var( yw (l)) = ivar(z i( l)). (A9) 



in probability as T ->• oo. By similar arguments, we have By Slutsky's theorem [74] and with (A6), (A7), (A8), and 

M 



that y E*=i ViMT^P) ~^ Var(y ii M,r(l)) = igVar (z,(l)) (A9), we finally obtain that 



ELl ^,M,t(%,M,t(*) -VtU ELi W<,M,t(*)) Ef=l Vi,M,T{t)) 

-—covv{x iM ,t,x 3M ,t) = 1 = (A10) 

V yg( M )T T,t=l(yi,M,T{t) - Vi,M,T) 2 T Et=lW,Ar,r(*) - Vj,M,T) 2 



converges in distribution to a standard normal random in probability with T e g(M) = ^pjy- 
variable as T — > oo. This completes the proof. 



Lemma 2 

For T — > oo, i? — > oo 



, M T -> 2$(-0) (All) Proof. With Lemma 1, we have that 
y/T e s(M) I 
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E 



H 



(r) 



ij.ALT 



y/T oS (M) 



= P Pij,M,T > 



P COrr(Xi,M,T,Xj,M,T) > 



9 N 
y/T cS (M) / 

= P 



+ P COrr(Xi,M,T,Xj,M,T) < 



T 



g(M) 



Pij,M,T >0]+P 



T 



g(M) 



Pij.M.T < -0 -> 2$(- 



as T — > oo. Furthermore, H\ r -' MT is bounded by and equality. 

1, so Var (^H^\j < J. By the independence of the R 
random networks 



(A12) 

as R — > oo. The lemma follows with the Chebyshev in- 
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